; docformat = 'rst'
;
;+
;
; :Purpose:
;
; :Inputs:
;
; :Requires:
;
; :Keywords:
;
; :Outputs:
;
; :Example::
;
; :History:
; 	Written by: Matt Rigby, University of Bristol, Apr 7, 2013
;
;-
function a2i_emissions, parameters, measurements, emissions_bias=emissions_bias, $
  emissions_no_modify=emissions_no_modify
  
  compile_opt idl2, hidden
  
;  Catch, theError
;  IF theError NE 0 THEN BEGIN
;    Catch, /Cancel
;    bad_io:
;    Help, /Last_Message, Output=theErrorMessage
;    print, 'ERROR READING EMISSIONS FILE: ', parameters['POLLUTANTS', pi]
;    FOR j=0,N_Elements(theErrorMessage)-1 DO BEGIN
;      Print, theErrorMessage[j]
;    ENDFOR
;    message, 'AIF_EMISSIONS: check emissions files'
;  ENDIF
;  on_ioerror, bad_io
  
  emissions=hash()
  Q=hash()
  Q_error=hash()
  Q_growth_error=hash()


  for pi=0, parameters['NPOLLUTANTS']-1 do begin

    nYears_species=(parameters['ENDY'] - parameters['STARTY']) - measurements['START', pi]/12
    q_Pollutant=fltarr(4, 12L*nYears_species)


    ;EMISSIONS ERROR FUNCTION, IF REQUIRED
    if keyword_set(emissions_bias) then begin
      nEmissions=12L*nYears_species
      P=fltarr(nEmissions, nEmissions)
      for i=0, nEmissions-1 do begin
      for j=0, nEmissions-1 do begin
        ;10 year autocorrelated error
        P[i, j]=exp(-1.*sqrt((float(i) - float(j))^2)/(10.*12.))
      endfor
      endfor
      
      seed=systime(/seconds)
      L=double(P)
      choldc, l, ld, /double
      for n=0, n_elements(l[0, *])-1 do l[n:-1, n]=0.
      l+=diag_matrix(ld)
      scaling=L##randomn(seed, nEmissions)*emissions_bias[pi] + 1.
      
    endif else begin
      scaling=replicate(1., 12L*nYears_species)
    endelse
   

    ;READ EMISSIONS FILE
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    record_start=0
    openr, fun, a2i_filestr(/input, parameters['CASE_NAME'] + '/emissions/' + parameters['POLLUTANTS', pi] + '.csv'), /get_lun
      comment=1
      str=''
      while comment do begin
        readf, fun, str
        if strmid(str, 0, 1) eq '%' or $
          strmid(str, 0, 2) eq '"%' or $
          strmid(str, 0, 1) eq ';' or $
          strmid(str, 0, 2) eq '";' or $
          strUpCase(strmid(str, 0, 4)) eq 'YEAR' or $
          strUpCase(strmid(str, 0, 5)) eq '"YEAR' then begin
          record_start++
        endif else begin
          comment=0
        endelse
      endwhile
    close, fun
    free_lun, fun

    data=read_csv(a2i_filestr(/input, parameters['CASE_NAME'] + '/emissions/' + parameters['POLLUTANTS', pi] + '.csv'), $
      header=header, record_start=record_start)
    q_year=data.(0)
    q_total=data.(1)

    q_semi=transpose([[data.(2)], [data.(3)], [data.(4)], [data.(5)]])

    if total(q_total - reduce(q_semi, /total))/total(q_total) gt 0.01 then stop
    
    ;EMISSIONS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    for bi=0, 3 do begin
      q_pollutant[bi, *]=interpol(q_semi[bi, *], $
        q_year, findgen(12L*nYears_species)/12. + parameters['STARTY'] + measurements['START', pi]/12., $
        quadratic=0)
    endfor
    
    if ~keyword_Set(emissions_no_modify) then begin

      ;MAKE ALL EMISSIONS NON-ZERO
      ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
      for bi=1, 3 do begin
        if mean(q_pollutant[bi, *], /nan) eq 0. then begin
          q_pollutant[bi, *]=0.01*q_pollutant[bi-1, *]
        endif
      endfor
  
      for bi=0, 3 do begin
        wh=where(q_pollutant[bi, *] lt 0.001*max(q_pollutant[bi, *], /nan), count)
        if count gt 0 then q_pollutant[bi, wh]=0.001*max(q_pollutant[bi, *], /nan)
      endfor
      
    endif

    ;Store emissions
    Q[pi]=q_Pollutant

    
    ;EMISSIONS ERROR
    ; If semi-hemispheric errors present in csv file, use that, 
    ; else use error from parameters file
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    if n_tags(data) eq 14 then begin
      q_Pollutant_error=fltarr(4, 12L*nYears_species)
      q_semi_error=transpose([[data.(10)], [data.(11)], [data.(12)], [data.(13)]])
      for bi=0, 3 do begin
        q_pollutant_error[bi, *]=interpol(q_semi_error[bi, *], q_year, $
            findgen(12L*nYears_species)/12. + parameters['STARTY'] + measurements['START', pi]/12.)
      endfor
;      for ti=0, 12L*nYears_species-1 do begin
;        ii=interp_indices(q_year, $
;          float(ti)/12. + parameters['STARTY'] + measurements['START', pi]/12.)
;        for bi=0, 3 do q_pollutant_Error[bi, ti]=interpolate(reform(q_semi_error[bi, *]), ii)
;      endfor
      Q_error[pi]=q_Pollutant_error
    endif else begin
   
      ;Absolute error
      q_error_pollutant=q_Pollutant*parameters['EMISSIONS_ERROR', pi]
      if parameters['EMISSIONS_MINIMUM_ERROR', pi] gt 0. then begin
        cutoff=total(q_error_pollutant[*, *])*parameters['EMISSIONS_MINIMUM_ERROR', pi]/$
          nYears_species/12.
        for bi=0, 3 do begin
          wh=where(q_error_pollutant[bi, *] lt cutoff, count)
          if count gt 0 then begin
            q_error_pollutant[bi, wh]=cutoff
          endif
        endfor
      endif
      Q_error[pi]=q_error_pollutant

      ;Growth error
      q_growth_error_pollutant=q_Pollutant*parameters['EMISSIONS_GROWTH_ERROR', pi]
      if parameters['EMISSIONS_MINIMUM_GROWTH_ERROR', pi] gt 0. then begin
        
        for bi=0, 3 do begin
          cutoff=total(q_pollutant[bi, *])*parameters['EMISSIONS_MINIMUM_GROWTH_ERROR', pi]/$
            nYears_species/12.
          case parameters['FUNCTION_Q', pi] of
            'X': begin
              q_growth_error_pollutant[bi, *]=sqrt((q_Pollutant[bi, *]*parameters['EMISSIONS_GROWTH_ERROR', pi])^2 + $
                replicate(cutoff, n_elements(q_pollutant[bi, *]))^2)              
            end
            'LN(X)': begin
              q_growth_error_pollutant[bi, *]=sqrt((q_Pollutant[bi, *]*parameters['EMISSIONS_GROWTH_ERROR', pi])^2 + $
                replicate(cutoff, n_elements(q_pollutant[bi, *]))^2)/2.7
            end
          endcase          
        endfor
        
      endif

      q_growth_error[pi]=q_growth_error_pollutant/q_pollutant

    endelse
    
  endfor


  emissions['EMISSIONS']=Q
  emissions['EMISSIONS_ERROR']=Q_ERROR
  emissions['EMISSIONS_GROWTH_ERROR']=Q_GROWTH_ERROR

  return, emissions
  
end